#########################以下为针对每个样本算ASE
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan")
sel1=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=c("X2B_X1T","M8_M7","M6_M5")
sel1=c("M2_M1","M48_M47","M50_M49")
sel1=c("M18_M17","M20_M19","M22_M21","M40_M39")
tag=1
sel1=paste0(sel1,".snp")
library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)

bayes_p=function(ref,var){
    df=data.frame(y=var,Ntrials=ref+var)
    null <- logit(median(df[,1]/df[,2]))
    formula = y ~ 1
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
    m=m1$marginals.fixed[[1]]
	coef1 <- m1$summary.fixed
	beta0=coef1$mean
	lower_beta0=coef1$`0.025quant`
  upper_beta0=coef1$`0.975quant`
    lower_p1=inla.pmarginal(null, m)#lower_p
    upper_p1=1 - inla.pmarginal(null, m)#upper_p
    df=data.frame(y=ref,Ntrials=ref+var)
    formula = y ~ 1
    m2 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
    ms=m2$marginals.fixed[[1]]
    lower_p2=inla.pmarginal(null, ms)#lower_p
    upper_p2=1 - inla.pmarginal(null, ms)#upper_p
	rtmp=paste0(2 * (min(lower_p1, upper_p1,lower_p2, upper_p2)),"_",beta0,"_",lower_beta0,"_",upper_beta0)
	return(rtmp)
}
for (i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t") 
file=file[file$exits=="T",]
result=matrix(,nrow(file),ncol=8)
if(dim(file)[1]!=0){
fn="normal_tumor"
for(k in 1:2){
files=file[,c(paste(unlist(strsplit(fn,"_"))[k],"_reads1",sep=""),paste(unlist(strsplit(fn,"_"))[k],"_reads2",sep=""))]
for(j in 1:nrow(files)){
reads=as.numeric(as.matrix(files)[j,])
ref=reads[1]
var=reads[2]
tmprt=bayes_p(ref,var)
	result[j,4*k-3]=as.numeric(unlist(strsplit(tmprt,"_"))[1])
	result[j,4*k-2]=as.numeric(unlist(strsplit(tmprt,"_"))[2])
	result[j,4*k-1]=as.numeric(unlist(strsplit(tmprt,"_"))[3])
	result[j,4*k]=as.numeric(unlist(strsplit(tmprt,"_"))[4])
	}
}

file$normal_bayes_pvalue=result[,1]
file$normal_bayes_beta0=result[,2]
file$normal_lower_beta0=result[,3]
file$normal_upper_beta0=result[,4]
file$tumor_bayes_pvalue=result[,5]
file$tumor_bayes_beta0=result[,6]
file$tumor_lower_beta0=result[,7]
file$tumor_upper_beta0=result[,8]
file_name=paste0("../bayes_pvalue_beta0/",gsub(".snp","",sel1[i]),".bayes_p.txt")
write.table(file,file_name,quote=F,row.names=F,sep="\t")
}
}